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SUMMARY 


An investigation of approximate theoretical techniques for predicting aerodynamic 
characteristics and surface pressures for relatively slender vehicles at 
moderate hypersonic speeds was performed. Emphasis was placed on approaches 
that would be responsive to preliminary configuration design level of effort. 
Supersonic second order potential theory was examined in detail to meet this 
objective. Shock layer integral techniques were considered as an alternative 
means of predicting gross aerodynamic characteristics. 

Several numerical pilot codes were developed for sample three dimensional 
geometries to evaluate the capability of the approximate equations of motion 
considered. Results from the second order computations indicated good agreement 
with higher order solutions and experimental results for a variety of wing 
like shapes and values of the hypersonic similarity parameter M<S approaching 
one. 
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1. INTRODUCTION 


Progress over the past 20 years in hypersonic aerodynamic technology has 
emphasized space re-entry vehicles. Relatively little advancement has occurred 
during the same time period in similar applications to hypersonic cruise aircraft. 
At lower speeds, significant gains in supersonic cruise efficiency have been 
realized using linear theory optimization techniques. Similar improvements appear 
feasible for slender hypersonic vehicles. Recent studies have shown that, at a 
Mach number of 6.0, both passively and actively cooled airframe structures with 
relatively small leading edge radius can be advantageous ly employed in hypersonic 
aircraft optimization. To characterize the aerodynamic implications of such nose 
shapes and evaluate the consequences of proposed beneficial interference concepts 
such as those involving propulsive-airframe integration, new analytical tools 
are required. Although substantial capabilities exist using Euler solution in 
conjunction with floating shock fitting methods, quicker response analysis tools 
are required to derive and optimize base point configurations. Less exact non- 
linear theoretical formulations hold the promise of meeting this objective and 
providing economic design codes which are responsive to preliminary vehicle 
definition efforts. 
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2. LIST OF SYMBOLS* 


c 

mean aerodynamic chord 

CD 

drag coefficient, D 

q s 

CL 

lift coefficient, L 

q s 

Cm 

pitching moment coefficient 

Cp 

pressure coefficient, P-P OT 

qoo 

D 

drag 

L 

lift 

M 

Mach number or moment 

P 

static pressure 

q 

dynamic pressure 

s 

reference area 

a 

angle of attack 

5 

flow deflection angle 
SUBSCRIPTS 

L 

due to lift 

00 

free stream 


*Additional specialized nomenclature is defined in the theoretical sections 
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3.0 METHODOLOGY 


Emphasis has been placed on approximate theoretical approaches which are 
capable of treating relatively general three dimensional problems but still 
sufficiently simple to be responsive to vehicle preliminary design efforts. 

The basic intent of the methodology is to produce future improvements in lift- 
drag ratio of hypersonic cruise vehicles. As a result of the strong impact 
that favorable interference has had on supersonic design and the use of such 
concepts in recent advanced hypersonic aircraft studies, candidate analysis 
should be general enough to systematically treat such problems. Finally, 
interest in high aerodynamic efficiency usually emphasizes relatively slender 
configurations at modest angle of attack; that is moderate values of the 
hypersonic similarity parameter. 

Potential theory was selected as a candidate methodology on the basis 
of its success to markedly improve aircraft efficiency at supersonic speeds, 
the ability of linear theory to predict the aerodynamic characteristics of 
slender vehicles at moderate hypersonic speeds (figure 1) and its successful use 
by Tsien in deriving the hypersonic similarity rule. The systematic extension 
of first order theory proposed by Van Dyke (1) was specifically chosen for 
further study on the basis of its ability to approximate exact two dimensional 
and conical flow results. See figures 2 and 3. 

Hypersonic small disturbance theory was selected as a second candidate 
methodology in recognition of the progressive non-isentropic behavior of the 
flow as the value of the hypersonic similarity parameter increases. Finite 
difference analysis of this approximation by Gunness (2) indicated that the 
solution was essentially as complex as that for the full Euler equations and 


3 






















EXACT 

FIRST ORDER 

SECOND ORDER 







thus would not be particularly responsive to preliminary design level of 
effort. Conservation law forms of the equation of motion applied over finite 
volumes will be examined as an alternative means of developing a solution. 


Integral theorems 
evaluated for the 
variations of the 


for predicting the forces and moments will be derived and 
rapid evaluation of gross aerodynamic characteristics. Simple 
flow properties through the shock layer will be utilized. 
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4. SECOND ORDER POTENTIAL FOPJ ULATIGN 


The equations for the invlscid flow of a perfect fluid are: 


continuity VC/ 3 ^) * o 


momentum 


energy 


- 7( w-’S.) 


C P T • 7 S..& 


- — Vi? i- It x (yxit) 

r 

c p r - ' t( u* - vO 


assume that at ©o 


LA. 


U*. £* * vn/„ e a 


and since 


V-l 


a. •* VRT 


the energy equation becomes 


t t 

Ov. = Ok. _ 


V-l 

z 


W„) -<£• It)J 


if 


V x L*. « Vs 


since 


o.*,(— ) 


the momentum equation becomes 


-7- V(u- IX) * ~ ~* Vf 5 


The continuity equation can be written in the foim 

t ex 

o- Y - w. -*■ — • u. ■> o 

or substituting from the momentum and energy equations 

~ W^)-(X w.)] j ( v-u.) = 7 V(lX CL) • IX 
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Since we have assumed yxl£.*o we can write (see firure 4) 


+ * 


S ■- w „ * 21 


<■ 


<5* , l£. X . jr 

3 * » a ‘ 7T 


{ a-1 ■ - ^ [ 2 < f, U.* #, W_ ) * i \ $ \ §( ] j ( ^ * s 


- i (u + { )-— - * <& — -»- 

L “ * ' 3X A 3^ 


C w -* £, >^r j j ^ u„* £ 2 w„* -j ( 2 } £^+ £) J 


Now we as stone small perturbations and introduce a small parameter £ , as 
well as an angle of attack £<*•» 




o„.+ w; <? = u„. y i 


2 l 

O'- M- 


U 1 * W 1 

«• #» 




t V ** ) * ^ t+ <V *» ]] C <U + ^ $ tt ) 
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Figure 4. Coordinate System Notation 
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or 


/T^vl [i ♦ a + r i- ~> eM r _.i <0 l . 2n 1 ^ 

V * L *■ <l+ocV}-* Xx i:3k L 21 J --/iTkV^ax 

= «i j [ 2 C <P X + «te <?, ) - 4> x * %+ ] [ ( l- rti ) ^ ^ ] 

* * 0<>a +c ?i] t, 

+ Kf* 

+ + 4,)^ ](«?*♦ <+ <9* ) j- 




«o u) 


now let <? * € ^ 

with the corresponding perturbation velocities 

Substituting into the above equation and equating the coefficients 
of each power of e we obtain the equations for the perturbation potentials. 

The first order equation is: 

0-0 <c* + M - r . o 

** 31 84 

The second order equation is: 


, <» 


(l - * 4> { 

m *» 31 r t ft 




For a planar problem the solution may be broken into lifting and 
thickness solutions for each order. Denoting a subscript <r for the 
thickness and p for the lift the equations become: 

first order 


CO 


<1-/0 * <p 0) * 6 

^ \* 


C I - Pj ) 4 >^ * * 9 ^ ' +• 
second order 


T“ * B* ’S' 

“ 11 *4 


f +• <p + <k 

p p V 

»» 31 »* 


/i M l n „c u> ,a 

? + 0+56 

<r <r r 

»K 


31 


r 


= n - £ {t 1 * V *?> *<£* > ’ ti- ((C’-j)} 

* * J 1 * * J 
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BOUNDARY CONDITIONS. 


If is the local normal to the surface of an object placed in the flow 
then the boundary condition is: 

{ ^ a ~ * ° 
where u. is the dimensionless perturbation velocity^ u. « ( , 9 a , «f ) 

If the surface is of the form 2 » 


then * /' 


€ s 


€ 3X c « “ * 3^ * e a 


and therefore 


£ | ~^==r * ^CXj^cZ<%a>] 


[*^,*^*,*>2 * - J ^= 

V« T < 


Now introduce the perturbation velocity expansion 


«? * c 4> 0i ♦ e * <j> a ' * 


v 

iO 


*i * 


_ oi , c») 

€ U. ♦ € U. f 


tM i 

< V * € V 


CM t U) 

c u: * c uj * 


and expand from the Z 3 0 plane 


U -[* J € Zcx.ai] 


^[<,3, ^ £o,a)] 


€ u^x^a, * € t u^(* J a,c£oc J 3)} * • 

c *^ 0 C*, a, o) * ^W/O] ♦ • 

Cl) if ci) . 4_ a fo , 

C u. 4*, 3,0) * € l*. <x t Z,°y- 4C<,a)^-W CX,3^0)]+ . 

€ W tl '(«,3 J 0) * € t [u3 U 'c< > i,®)*-ZC < »a)^‘-3’c* J 3 J o)J* 
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Equating powers of £ we obtain the boundary conditions for the various 
orders of the perturbation velocities at the Z = 0 plane. 


£•> , 

us C -a J o ) 


- «<„ 


o) 


<* w Cn o> . O) 

^cx.a) y cx.-a^i 


•_ , 3 to 

- ZC*.*> — ^ Cx,^ j0 ) 


From the first order equation. 


a «> 


" ~il w 3 0-0-^- ulV^o) + ~ v'u.^o) 


Then assuming an upper and lower surface 


2 * o 


Z = Cc*,a) + Tt-*,'. a) - «<.C.x- x„) 


2 * o 


2^.a) 


Ccx,^) - Tu,a) - * 0 (<- x #- ) 


the following equations result. 


| (I) + 0} -* ^ 

Y L w ^,i,o ^ ^v')j = 


T [ u>C ^x,-a j o*) - ui^'cxj^o*)] = TUja) 


uAlx^o')] = 


— c«,a; - C'Ax.-s,o') s 


+ -^(Vcx.^y ) ? V^Cx^O] 


♦ 0-M*) ^ -£L u - C ) t’S3,o t ) i ^’ ) Cx 1 ^ o-j] 


Li-nl) Tc*ji) — 


+ [Cc*,a>- «< 0 (x-x.^] — -r-A'o^y) * 'A^oj 


S t r U) + U1 v -t 

Tc* — tL ^ <•*,*, 0 ) ; V C^hjjoOj 

A -4 
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and since the problem is planar 


7 L U)tt 'cx,a,o*> ♦ 

a> •, 

° C*» a, o ) J « 

* a> 


«0*V. *, o') ] . 

^ C*. 3,0) 

* 

w ' C, 't X j ] 5 

0) 


uf'cx, 3 , 0*3 3 - 

t *3 

«> C^.a.o) 

i[v°W.«0 * 

v'V,v'J] - 


-jt V C,> Cx, 3j o*) - 

y^cx,^ o') ] = 

v> x <V) 


where P and <r refer to velocities due to lifting and source elements. 
In this form the boundary conditions become: 


<0 




<• X, a, o ) 
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The thickness distribution can be written in the following form. 


where 


where 

line. 




cca) t ( i J 3) 


a 


x- x te c^) 
local chord 


or fraction of chord 


* maximum value of t/c 


T^ c; ^ 


- < T u.r t t „) 






T 


T is the tangent of the sweep of the local constant percent chord 
Therefore. 


3 * 


TCx.,tO 


<=<»!> 


tc^ —■ 


- T i- 


l rca) 
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PRESSURE COEFFICIENT 


The pressure coefficient is given by the equation 






Since to second order the flow is isentropic 


JL. 

*-» 


Cp = 
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7^ * 


therefore to second order 
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These velocities must be evaluated on the surface 


E » c = Cc*,*) - Tcx jS ) - * 0 (x-x. 0 ) 


€ U.O., a, * -- 


Therefore to second order 


to x r w a f, » . ci> l 
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+ 2 uo ( °c Xj-ijO) | 
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Three different approaches were examined for solving the supersonic 
second order potential equation of motion. They were 


a) Discovery of a particular integral for three 
dimensional flow 

b) Source volume analysis 

c) Finite difference analysis 


4.1. Particular Integral Existence 


Second order theory expands the exact potential <5 in terms of thickness 
( S ) or angle of attack ' '• 

£ * U £ x * s 4> oi + & } 

j> u \ satisfies the usual linearized theory, </ M is the second order correction 
and satisfies a forced equation 


i"' 


q ¥ 


•C * C - fit 




/ 3 > * n„-i 


♦ 

aa 


~ r »* 


1 3 ff. Y-l l 7 .W 1 ,<’) Z 1 r 


An algebraic integral would have the form 






O * = 0 F 


CD 


Van Dyke has found such integrals for 2-D and axisymmetry. However, a 
straightforward study of conditions on F in the general case shows that 
no such algebraic integral can exist. To do this consider <£°= f given; then 
for example 





F f 



Comparing like terms on the RHS and LHS of Cl) gives a set of relations 
with no solution F. 
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4.2. Source Volume Analysis-- - 

Fcr a thick lifting wing system, the first order 
solution may be obtained using a Woodward type finite element. Using this 
scheme the entire planform is divided into quadrilateral panels having two 
streamwise edges and placed in the mean chord plane. Each panel has a 
distribution of source and vorticity strength to account for the effects of 
thickness and lift. The source and vorticity singularity strengths are 
, determined by satisfying the boundary conditions at a set of control points. 
When these strengths are known the first order velocities may be computed 
anywhere in the field. In particular the value of 


( c* 


used in the solution of the second order velocity potential can be calculated 
anywhere. This function is computed at the centroid of a set of source 
volumes distributed throughout space, giving the strengths of these source 
volumes . 



Using the influence equations for these source volumes , the 

velocities induced back on the planform from the spatial source distribution 

may be calculated. The source volumes have the property that: 

P within volume 

£ = p is a constant 

o outside volume 
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The source volumes to be considered have top and bottom faces 
composed of two identical quadrilateral panels separated by a distance 42. 
The panels have two streamwise edges and swept leading and trailing edges. 
The corresponding comers of the panels lie along lines x-n = const. 
Therefore the two side faces have two streamwise edges and leading and 
trailing edge sweeps of magnitude Z. 



£ * § * £ 
** 33 


inside the volume 


outside the volume 


where 


§ x + z - * <$> 


The 4>i refer to the velocity potential contributions from each of the 
four panels which form the top, bottom and sides of the volume, and <p is 
an additional term which is nonzero only within the volume formed by 
the four planes of the top, bottom and side faces. 
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The 4* for each of the four panel faces is composed of a contribution 
from each comer, each corresponding to an integration limit. 



, On any given panel or face, at a comer with an edge having a sweep T and 
where 'u is the sweep of the edge at the corresponding comer on the 
adjacent panel 


4> 


where 


It* + *'*/*'} 


l 4+ * 4,1 


/b - | - ri* A.* 

rN 


ni<l 




u. = — 

* ax 


-- 




& - 
y 9i 


~^r | T±i t - (tV)z^ + (x-T^)-f 3 J 
_L ♦ TKrVJf, - TCx-r^f,- i } 

{ [TC^-Tai- a 3 V - (rV)?i 3 ^T(^f T }^ 


“’••-ST' s - + j 

v - - 4f- 1 - - TU * T ^ - - * j 

Hs s -§r~ : ~ lk { Tz *. " 2 K + j 

1 , 1 n R »* r . 1 i n R ♦ CT»4/3a) ^ gR 

1 2 R ' X 1 /F\T* 1 >i4r ^V T ‘ > /3‘ !*. ‘-3 ' 


The velocity potentials 4>^ and have the property that 


ax 




where 4>v m and Aa o are the velocity potentials induced by panel comers on 
panels having constant vorticity, and constant source strength respectively. 


22 



The additional term <p is such that 


outside V t 


= O 


inside 


*4 


1 

AO 

a x 

t-rSrV'] 2 l 

34 

- <r A 1 

! T, 11,1 - T 
■ 

*3 

‘ tr 4 + r*r/3 , i 2 | 

a«P 

- <r- -8- 

[ is.i - rjs 4 

32 

Lt** z 


where V* is the volume determined by the planes comprising the four panels, 
and -wv. . t 

S, - - t. - 

i*. * *•» * *7*v] - r i la- *[*-iu,*s»)] 

♦ 

where and are the leading and trailing edge sweeps of panel I. Therefore 
and i 2 are the x distance behind a point on the front and rear faces of 
the source volume which has the same y and z values as the point being 
influenced. For supersonic flow 5 ^ and 5 2 are zero upstream of their 
respective faces. 
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Each comer yields a contribution to the volume influence coefficient from 
each of two adjacent faces. This fact may be used to write the influence 
equations for a given comer in a different form. On a panel whose normal is 
in the z direction the influence coefficients yield. 


’ u.' 

• V > 

UJ 


F (*» ■fcjT, T ) 


On the adjacent face, with normal in the y direction, the corresponding influence 
is 


’ IA. 
• UJ 

-v 


- F ( x , Z,^,T,-T) 


where (x,y,z,T,r) and (u,v,w) are both expressed in the coordinate system of 
the face with the normal in the z direction. Therefore, the equivalent set 
of influence coefficients for a comer of a given face are: 


u. = 

I 

zrre ] 

L - 


V = 

<rd«. 

2ire 

^ ~ €Z-f, Ta«-^ 4 


u? = 

J 

Zi re 

L ' 



€ 

-- t' + p 

i - *-T^ 


1 

} 


- 


- ra 


The (x,y,z,t,T) and (u,v,w) are now expressed in a coordinate system whose origin 
is at the comer and whose z axis is in the direction of the normal to the face. 
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As € * T^-r'vs*— o the quantity ( ‘r-f t - -f^) will cancel with a corresponding 
term on the adjacent face and the velocities become indeterminate. By expanding 
about € = 0 the following form for the velocity influence coefficients results. 



This form is indeterminate when T = 0 and must be rewritten in the following 

form. 

T s O 6*0 


IA. = 


V = 


u> 



- af . * T * F , -F } 

♦ + T [ T »^- RH-} 

- fih - } 


m - o 6*o 


IA. = 


V = 


u) = 



<rA 

2.7T 


l 


— 2 "f + — - Z — — 1 
2 2 jr j 

■ i 5 h ■ } 
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'Jliis spatial source distribution can then be used to solve for <^(2) . 


a <f> - + a 

A solution for <£^can be obtained by using a solution to the equation 


t C2> rw \ 

aa ea 1 



o 


with the boundary conditions adjusted to cancel the normal velocities induced 

on the planform by the spatial source distribution, f> (x,y ,z) . Therefore 

4 ( 2 ; 

on the planform may be obtained by using the same Woodward panel scheme 
with boundary conditions determined by the normal velocities induced by 
p(x,y,z) and the second order boundary conditions. 


A pilot code based on this approach was developed. Computations and 
comparison with higher order solutions and experiments is presented in 
section 6.1. 
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4.3 Finite Difference Analysis 


The first and second order theory equations and boundary conditions derived 
in the previous sections are rewritten here for completeness with a slightly 
different notation. 


First Order 


□<£, = + <}>, - 0 2 q, = 0 

yy zz xx 

with boundary conditions on the body given by 

3F 

<f 1 (x,0,z) = (x,z/b) - a 


(1) 


( 2 ) 


Second Order 


U H ~ *2 + *2 


B 2 <t> 


= M 2 4 fl A? + d> 2 


yy z: 
with boundary conditions 


2 xx 3x 


2 ^ n x n y 


4 f > 2 (x, 0 ,z) = (x,z/b) - aj (x, 0 ,z) 

- |f(x,z/V) - ax j 4 ^ (x,0,z) 


(4) 


+ b' 1 <*>, (x,0,z) (x,z/b) 


yy 

3F 

sW 


where y = 6F(x,z/b) + ax defines the body geometry, and a = a/6. For a flat 
plate delta wing, a = -1 on the compression side, and a * 1 on the expansion 
side. The angle of attack is a. 

The pressure coefficient is computed from 

P-P / \ 

= - 26$ - 6 2 2q 7 - s 2 ® 2 + q 2 + b 2 . 

^°° -v \ "x x y z / 


(3) 


(3) 
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Finite Difference Procedure 

In this section two different finite difference procedures are discussed. 

In one, the velocity potential is determined from a single equation and in the 
other, the velocities u = 4, v = <£ v and w = are computed by treating a svstem 
of partial differential equations for these quantities, denoted hereinafter as 
the systems approach. 

Scalar Approach 

Equations (1) and (2) are hyperbolic m the x-direction provided 
S 2 = - 1 is positive. Considering a grid molecule as shown in figure 5, for 

an explicit finite difference procedure the information at point (i+l,j,k) can 
be computed based on past information. A simple explicit procedure for 
Equations (1) and (3) are obtained by central differencing all the terms m the 
equation about point ( l , j , k) . Thus, for equal mesh spacing (Ax, Ay and Ac 
are constants), the appropriate discretications are: 
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The term ^ J . is differenced similarly. All these finite differenced 

quantities are substituted into Equations (1) and (3) and the first order 

velocity potential <b, and the second order velocity potential 

i+l,j,k “i+l, j,k 

are solved for. The boundarv conditions are imposed using a boundary point 

operator. 


Boundarv Point Operator 


The boundary conditions can be applied m two different ways. Figure 6 
shows the two arrangements. In one (refer to figure 6a), grid points are placed 
right on the body and dummy points below the body are used. For this arrangement, 
the term at the body point (point 1) is differenced m the following way. 

$2 * 

5v* • 
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Ay 

0 9 -K 


iM^ire 6a. Dummy Point Arrangement 



Mgure 6b. No Body Point Arrangement 



The dummy point value bg is replaced m terms of the body boundary condition 


(b v ) as 


*0 = ^2 " ZAy 


at 1 


Substituting for bg we get 


CO 

y at 1 


2b, 


_ 9 


b, - 2Ay (bj 


T at 1 


(Ay) : 


In the other arrangement (refer to figure 6b) , grid 
the body. While differencing for b at the first point 
body boundary condition by will come in to play. 


CO 

77 at 1 


4> 7 - b n - (b v ) Ay 
= 1 7 o 

[Ay]^ 


points are not used on 
above the body, the 


Stability 

In general , all hyperbolic explicit schemes are subjected to stability 
conditions to achieve convergence. This stability condition is usually in 
terms of the mesh spacing Ax, Ay and Ac. Basically, stability is achieved if 
the discreticed domain of dependence completely includes the continuum domain 
of dependence. This is illustrated m figure 7. The continuum domain of 
dependence for the point P(i+l,j,k) is shown by the dotted line Mach cone whose 
projection on the i tn plane is represented by points B',C',D' and E’. To 
ensure stability, the grid points B(i, j ,k+l) , C(i,j-l,k), D(i,i,k-1) and E(i, j+1 ,k) 
should fall outside of Mach cone projection B’C’D'E'. This means that the discreticed 
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\tach rhombus BCDEBP must enclose the continuum Mach cone B'C'D'E'B'P. This is 
achieved if the following condition is met. 


1 / Ax 2 + Ax 2 \ 
fF \ Ay“ AF ) 


< 1 . 


For a special case Ay = At, the stability condition is 


Ax < _s_ 
Ay /j ‘ 


C6) 


Usually, best results are obtained when the equalitv sign m Eq. (6) is satisfied. 


Artificial Damping 


The explicit scheme that is used to solve Eq. (1) is completely neutral. 

If one solves for the spectral radius of the amplification term in the frequency 
domain we get 




y 


+ 



A 




r 


3/ 

) 2 1 

M 

1-2 

v 2 sm 2 

L y 

-f- + vi sin 2 

2 

!] 




For a wave solution, we require that the spectrum be complex, -which means the 
quantity m the radical must always be positive. This term is positive only if 




_ 1 Ax _ 1 Ax 

v y S Ay ’ v c 3 Az 


(3) 


which is the stability requirement on Courant number. However, the real reauiremen 
is that 



< 1 . 


(9) 
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We can see that the restriction (3) provides us with roots which always lie 


on the unit circle for all stable Courant numbers. Any discretization errors 
which are introduced can never damp out. They will simply be trapped in the 
mesh as the solution proceeds. Another interesting property of this scheme 
is that the usual damping terms which are appended to the right side of the 
difference equation are of no value. They are either destabilizing or fail 
to alter the modulus of the eigenvalues. 

To damp out any errors that creep into the solution, a smoothing operator 
to the final result has been introduced. Let be the difference operator 
such that 

*i+l,J,k = L C * 
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New lee a purely dissipative operator Lq be applied to the result and let 

this result be . such chat 

i+l.l 


final _ , 

i+l, j,k “ S ^i+1, j ,k ’ 


Then 


.final _ r t 3 
i+l, 1 ,k u C i, j ,k 

For stability, we require 

II L c V < 1 


( 10 ) 


( 11 ) 


( 12 ) 


or 


II LjjII II Ljl < 1 . 

Since IlL^II = 1 we simply require that IILqII ^ 1. Any operator wmch is 
dissipative will not disturb the stability of the system. It is suggested 
that the difference operator be of the form 


(1 -t 5 P *1+1,3, k 


(13) 


where the z-direction is indicated. The symbol o£ represents a fourth order 
differencing operator. A similar damping term in the y-direction is added 
with no problems. For this type of operator, the difference equation becomes 


.final _ cj 

i+l,3 A °i+l,3,k 3 


- L 


- 4o 


i- 1 -! 


i+l , J , k +2 1 - 1 , 3 ,k+l 

,3 + ^ 3 »’*“2 ] 


6o . , , 


( 1 *) 


Stability is assured since 


I Loll 


1 - 2 o sin 4 — j 


ir 


0 < u) < \ . (15) 

This form of dissipation should leave the low frequency components of 
the solution relatively unaltered while the hign freauency waves are 
attenuateu. The formal accuracy of tne solution snould net oe altered 
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using this procedure. The damping really appears as a fourth-order term and 
the second-order result obtained should still be intact. Solutions obtained 
with and without using this damping operator are discussed in Section 6.1. 

Systems Approach 

Although the potential formulation is fairly straightforward to solve with 
only one dependent variable, it has a few drawbacks which are described below. 

To alleviate some of these deficiencies, a different approach has been studied. 

Here, instead of solving the single velocity potential equation a system of 
equations in u, v, and w (perturbation velocities) are treated. 

Comments on Velocity Potential and Svstents Approach* 

1. The <f> equation in a Cartesian coordinate system involves second 
derivatives. This poses a problem if we have to make any coordinate 
transformation because the transformed equation may get so complicated that 
even a simple explicit scheme would require time and storage consuming 
matrix inversion procedures just as in the case of implicit schemes. 

2. To solve the problem in Cartesian coordinates requires too many 
field points. This can be avoided by means of a simple coordinate 
transformation (preferably body fitted system) . On the other hand, 

a velocity potential formulation is not very accessible to transformations 
of any kind for the above reason. 

3. The u,v,w systems approach involves only first derivatives, and can 
be put in conservative form for any arbitrary general coordinate 
transformation (described in the later section of this report) . 
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4. The existing three-dimensional Euler solvers use continuity, and the 

three momentum equations in vector form along with the integrated energy 
equation. The u,v,w systems approach has the strong advantage over the 
velocity potential formulation in the sense it can be easily incorporated 
into any one of the existing Euler codes. Within such a framework, the 
perturbation velocity formulation is linear and simple looking compared 
to Euler equations and hence can save computational time. This simplified 

described. 


( 16 ) 


(16) are the vorticitv equations. Here 

u l = <j>i , v^ = 0^ and * = where is the first order perturbation 
X V L 

velocitv potential. Eq. (16) can be written as 

E 1 " F 1 + G 1 = 0 * ( 17 ) 

x y t 

The first order boundary condition is 

v x (x,0,s) « F x - | (IS) 

where v = cF(x,c) + ox prescribes the body shape. The angle of attack is 
denoted by a. 


formulation will now be 
First Order Equations 



The second and third equations m 
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Second Order Equations 



M* . 
a 7- 3x 






u* + v* + w* 



= 0 


> (19) 



/ 


where u 2 = < J > 2 > v 2 = ^2 ^ w 2 = ^2 ’ ^ ^2 ^ seconc * orcler perturbation 

x y z 

velocity potential. 

Equation (19) can be rewritten in the following way : 


E 2 + F 2 + G 2 = 0 
x y z 


where 



F 


2 




6^ 

0 


u„ 



( 20 ) 


39 



The second order boundary condition is 


v 2 (x, 0,2) = (F x - f) Uj. - tF - f x) v x + Kj F. 


( 21 ) 


The first order equations (17) and (IS) and the second order equations (20) 

I 

and (21) are m Cartesian system. Any general transformation of the form 



( 22 ) 


does not change the form of these equations . After the transformation 
Eqs. (17) and (20) result in 


AAA 

E-. + F, + G, =0 

c n 5 


A A 


E~ + F- + G 7 =0 
C "n 


ivhere 


E = — 
C 1 J 


F l = 


E n n + F-, n + G. n 
l x 1 y 1 z 


G l = 


pr+F^+Gt 
E 1 '•x 1 ’v U 1 


(23) 

(24) 



E 9 n + F 7 n +G, n_ 
p = _r c “ — L r — - 


; E 2 S + F 2 V + G 2 

G 2 = 


and J is the Jacobian of the transformation given by J = q 5 y - n v 5 x - 

At this point, ive note the identical form of the first order Eqs. (17) and 
(23) and second order eauations (20) and (24). By contrast, m the velccitv 
potential scalar formulation any transformation completely alters the Cartesian 
form of the equation. 
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Conical Wing-Body Problem 


To study the usefulness of the systems approach, conical delta wings and 
wing-body combinations with supersonic leading edges are studied first. To 
take advantage of the conical nature of the flow field and to map the leading 
edge planar shock into a constant coordinate surface the following conical 
transformation is incorporated. Referring to the coordinate system of figure 8, 
we let 
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Corresponding to the transformation (10) we have 



and 



The computational domain (n,C) is a rectangle as shown in figure 9b. Various 
other transformations besides (25) are also possible. For example, ; = x, 
n = y/x, and E, = z/x is another such conical transformation. 


outational Boundaries 


M=l, L=l, L is the plane of symmetry of the delta wing. 

max 1 a 

M=M , L=l, L is the right boundary of the computational domain, 
max max = 

Physically, points on this boundary lie in the region where the solution is 
two dimensional. 

L=1 . M=1 , M is the body boundary or it is the boundary to which the 
boundary conditions are transferred and v«here the boundary conditions are 
applied. 

L = Lmax’ ‘^ I=1 ’ M max the t0 P boundary lying in the physical region 
undisturbed by the presence of the wing (cone of silence) . 
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MacCormack's Difference Scheme 


Applied to equation of the form 


AAA 

E + F + G,. = 0 

Z n c. 


pn+1 pn A? (On pn \ Ac /pn pn \ 

C L,M l ,M ' An \ L+1,M ‘ l,M/ ' AC \ L,M ' b L,M-l/ 


pn+1 _ 1 F n + F n+1 A ^ /pn+1 9 n+1 \ A ? /pn+1 pn+i \ 

L,M 2 C L,M ’ An \ L,M ' ' AC \a,M+l " L,M / * 


Once is computed (u,v,w)£ + ^ is known. For a conical problem as n gets 

large, the solution approaches a radially asymptotic value. The marching step 


size Ac is chosen to satisfy the stability conditions. 


Boundary Differencing of Boundary Conditions 


At the top boundary the values of u,v,w are kept fixed at zero (undisturbed 
flow region) . 

At the symmetry boundary the equations are integrated using a forward 
predictor and corrector for the differencing along varying M index and the 
boundary condition w=0 is applied. 

At the body boundary the equations are integrated using a forward predictor 
and corrector for the differencing along varying L index and the boundary 
condition v=prescribed is applied. 

At the right boundary which is m the region of two dimensional solution, 
it is assumed that 
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Computationally this is implemented by setting 


(i^, v r w r u 2 , v 2 , w 2 ) = (u 1( v r w r u 2 , v 2 , w 2 ) ^ 

%ax 1 'max 


Symmetric Subsonic Leading Edge 

In general, if the leading edge is subsonic, then the bottom and the top 
surface communicate which means the computational domain must include both the 
top and bottom right hand side quarter plane. The special case of zero angle 
of attack symmetric wedge delta wing case shown in figure 10a can be handled by 
considering only the right hand side upper quarter plane as shown in figure 10b. 
The boundary condition on the wing is just v=l and off the wing v=0. Calculations 
for such a case using a transformation ? = x, q = y/x, and £ = z/x are shown in 
Section 6.1. 
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a) Physical Domain 

Figure 10. Symmetric Delta Wing Subsonic 
Leading Edge Case 





5 . INTEGRAL FORMULATIONS 


In addition to the previous approach, conservation law forms of the equation 
of motion applied over finite volumes also appear promising for rapid turnaround 
aerodynamic analysis and design methods. Referring to figure 11, a control volume 
denoted as ABCD in the plan view will be considered. Without loss of generality, 
we will illustrate the method to be studied within the framework of 
hypersonic small disturbance theory' (HSDT) . Generalisation of these 
developments to the Euler equations is straightforward. The HSDT model should 
be adequate to describe flows for freestream Mach numbers, M greater than 5 
and characteristic flow deflections, A, generated by incidence angles cs, and 
thickness or fineness ratios 5 of the order of 0.2 radian, or more generally, 
flows with the hypersonic similarity parameter 1/M^A 2 of order unity. For 
higher values of H, or lower incidence combinations, a different type of integral 
method based on Prandtl Glauert and second order theory will be formulated and 
evaluated herein. 


5.1 Hypersonic Small Disturbance Theory 

The HSDT problem is obtained bv substituting into the exact equations 
asymptotic expansions which are approximate representations of the velocity' q, 
pressure P, and density' o in a limit involving a characteristic flow deflection 
parameter o, and M^. In this limit, strained coordinates are required as m 
boundarv laver theory to keep the flow field between the shock and bodv m view. 
The asymptotic expansions are 


§ Cx,v,2,M o8 ,5) = 


i (1 + 6 2 ufx,v, z,H) + ...) + 6v(x,v,z)i - Swk 


(la) 


48 



BOW SHOCK 



Wjjuro 11. Control Volume for Typical Vehicle Geometry 


(lb) 


P - P 

CO 


p U2 

C 


rs 2 

0 p + . . . 



(1c) 


for H = 1/M7$ 2 , x = x, y = 7/6, 2 = 7/5, B = b/6, A = a/6 fixed as 5 -*■ 0 inhere 
b is the maximum lateral dimension of the configuration, 5 is the fineness 
ratio, is the freestream velocity, °° subscripts refer to freestream conditions, 
i, j, and k are unit vectors in the x, y, and 7 directions, respectively. 

In what follows, attachment of the shock only at the nose will be assumed for 
convenience. Generalization of the analysis to more complicated situations 
will be considered in the future. 

Substituting Eqs. (1) into the exact equations of motion and retaining 
like order terms gives the HSDT equations 


a x + V T ' V T = 0 

(2a) 

+ V-pP = 0 

(2b) 

= 0 

(2c) 

where 


D 5 3/3x + V T • P T 


V_ = — j + — k 
T 3y J az 


Vj = vj -*■ wk 

(2d) 

'to 

a: 

in 


Y = specific heat ratio . 
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On substitution of (1) into the shock conditions, the following approximate 
relations are obtained on the shock S(x,y,z) = 0 


V T x VS = 0 (3a) 

|V T | = —j- CT - H/T] (3b) 

where 


T i S.,/17^1 (3c) 

p s 1 p| s»o = (TJ ■ H1 (3d) 

- pt * T ' 2 • w 

The body boundary conditions are obtained in a similar manner for flow tangency 
on the body B(x,y,z) = 0, giving 


B x + V T • v tS = ° . (4) 

Note that if S is known, (3a) and (3b) uniquely determine V^/S. 

One approach which has been considered utilizes conservation theorems and the 
HSDT formulation embodied m Eqs. (l)-(4). This method is the full nonlinear 
analogue of the integral theorems and area rules developed by Rockwell for 
hypersonic lifting wing body combinations in Refs. 3-5. There, a linear 
perturbation problem was obtained for configurations with "very supersonic" 
leading edges. Applying these theorems to the previously indicated control 
volume, we obtain 


A 


II adA 


A A 


, continuity 


(5a) 
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2L = . / 
p-u 2 s 4 

f crvdvdz 
/A 

, y momentum 

(5b) 

l 

/ owdydz 
i A 

, z momentum 

(5c) 

20 . - i 

T / p + H/y 

* § (V 2 * w 2 )) 

dydz + HAg/y , x momentum (5d) 


A U - 1 


where AA=A^ - Ag, L, Z, D, are, respectively, the lift, side force, and drag 
and Cg, C 7 , and are their normalized counterparts, and the area integrals 
are taken over the projection of the shock layer m the section A-A. Analogous 
expressions for the moments will be obtained using the moment of momentum 
theorem. 

To evaluate the integrals in Eqs. (5) various methods have been considered. 
One approach that could lead to a fairly rapid procedure is to make an "averaging" 
hypothesis that neglects variations m the flow quantities along lines such as 
OP in the cross flow plane in figure 11. Thus, the values of the integrands m 
(5) are taken along the shock. For example, if a polar coordinate system is 
used with the lines OP assumed as 9 = const., as shown m figure 11, then 
Eqs. (5) simplify to 


C L- 


/ Ka s v s d6 


. 2tt 


C 7 = - f KbcW-dG 

J Q Si 


^ - / 0 K 


p s + H/y o s (v| h- w |) HA b 
Y - 1 + - Y 


d9 


(6a) 

(6b) 

(6c) 
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where K = G 2 (x,0) - F 2 (x,0) with the explicit representation of the shock and 
body taken respectively as 


S = r - G(0,x) = 0 

, shock 

(7a) 

B = r - F(0,x) = 0 

, body 

(7b) 


where r 2 = y 2 + z 2 and 0 = tan -1 y/z. 

The integrands in (6) can be calculated from the shock relations (3) 
providing the shock shape S is known. In line with the previous assumption, S 
can be determined by equating the value of the transverse velocity on the body 
to its corresponding shock value and using the shock relation (3b) to obtain 
the following nonlinear first order partial differential equation for G: 

T - H/T = Xp Fiy VF 2 + F 2 (8a) 

where 

|V T I = FF^Vf 2 + F 2 (8b) 

s 

T = |V T | = GG^VG 2 + G 2 (8c) 

S 

and 

|V T | ■< VB = 0 (8d) 

1 S 

has been assumed. 

Equation (8) can be solved as an initial value problem using the characteristic 
construction described m Ref. 6, or other numerical methods. This procedure 
should be fast since there are only two independent variables involved. In 
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particular, Eq. (8) simplifies to an ordinary differential equation for axial 
symmetric or conical bodies such as delta wing. 

Note however, that Eqs. (8) and (3d) uniquely determine p s without the 
need to solve (8) as a differential equation for the shock shape G. If the 
additional approximation over and above those used is made that 

P s « P B ( 9 ) 


(consistent with thin shock layers at high Mach numbers) , then a surface 
pressure integration rather than the momentum theorems of Eqs. (6) can be used 
to determine the inviscid overall forces and moments. In analogy to Eqs. (5), 
these are given by 

c l ' pjn; ■ 2 // s % ^ tl0a: 

2 


c z - 5JFT - 2 //«. Ps "•* as 

7 

C D = FTF"s = 2 s n ‘ lds ^ 10c ^ 

=° s 

~1 


where S refers to the surface area of the configuration. 

For lower hypersonic Mach numbers of the order of 4 to 6, the assumption (9) 
is questionable, and the formulae (10) must be modified bv replacing p s by Pg, 
the true surface pressure. The relative accuracies associated with application 
of (5) as contrasted with (10) when certain simplifvmg assumptions are made 
concerning the structure of the shock laver or disturbance field are under 


54 



examination. If the conservation laws associated with (5) are to be used, then 
the factor K in the integrands must be determined and the differential equation (8) 
for the shock shape must be solved. 

Assumptions (8) and (9) appear to be reasonable approximations at high 
supersonic Mach numbers as shown in figure 12 taken from Ref. 7. Calculations 
for = 03 are indicated for a right circular cone at zero incidence where for 
lower Mach numbers these assumptions correspond to the wedge curve in the figure. 

A partial evaluation of the method was presented m Ref. 8 for axisymmetric 
cones at incidence and zero angle of attack. In that assessment, the ambient 
pressure was neglected compared to the static pressure on the shock layer. This 
provided agreement of the order 16-20% with the numerical solution of the Euler 
equations from Babenko et al. given in Ref. 9. A somewhat more consistent 
procedure which has been recently investigated and accounts for the ambient 
pressure increases the discrepancies to approximately 36% as indicated in 
Section 6.2. Similar inaccuracies will be discussed there in connection with 
studies on elliptic cones performed in this aspect of the contract. To reduce 
these errors, additional effort is proposed to construct some sort of nonlinear 
feedback model involving the equations of motion. This concept will be discussed 
again in Section 6.2. 

5.2 Prandtl Glauert and Second Order Theory — Transverse Integral Method 

To simplify the calculation of overall loads and moments and bypass some of 
the difficulties associated with pressure calculations according to second order 
theory, an integral method has been developed. Although similar m the objective of 
reducing the dimensionality of the problem, the thrust of this technique differs 
from that embodied m the integral methods described m Section 5.1. Rather 
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0 20 * O os- 


Figure 12. Surface Pressure for a Wedge and Cone 
According to the Exact Theorv for 
M = 00 and y = 1.4, and According to 
the Newtonian Model. 
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than dealing with volume fluxes of conserved quantities such as mass, momentum, 
and energy as Eqs. (5), we consider integrals of the pressure that specialize 
to loading along strips on the body approximately normal to the freestream 
direction. This "transverse integral" technique, originally developed by Lagerstrom 
and Miles, has been used by Malmuth in Refs. 3-5 to estimate L/D performance of 
conical wing-bodies at hypersonic speeds. In this section, we foimulate similar 
rules for Prandtl Glauert theory and outline the procedure to be further developed 
to treat the second order linear problem. 

With the notation associated with the winglike ogive shown in figure 13, we 
consider initially the transverse integral technique appropriate to Prandtl 
Glauert theory. Denoting u as the perturbation velocity in the freestream 
direction the appropriate equation of motion for the perturbation velocity 
potential <f> and boundary conditions when differentiated with respect to x the 
streamwise coordinates are 

□u = 0 (11a) 

Uy(x,0,z) = f^ 

where □ = 3 2 /3y 2 + 3 2 /3z 2 - 6 2 3 2 /3x 2 , 0 2 = - 1, and the body is given by 

the equation B = y - 6f(x,z) =0, with 5 a characteristic flow deflection angle 
on thickness. .Assuming for initial simplicity that the ogive has supersonic 
leading edges here and in what follows, a characteristic surface emanates from 
the leading edges. This surface is shown schematically in phantom m figure 13. 

We remark here that if the edges are subsonic, certain additional relations over 
those to be discussed must be introduced to treat the initially unknown velocitv 
components off the wing. 
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CHARACTERISTIC ENVELOPE 
(SUPERSONIC LEADING EDGES) 


Figure 15. Notation for Transverse Integral Method m 
Prandtl Glauert and Second Order Van Dyke 
Theories 
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Returning to the supersonic leading edge case, we define a spanwise 


"pressure" integral of the form 
z c (x,y) 

Q(x,y) = / u(x,y,z) dz (12) 

0 


where z = z c (x,y) is the equation of the characteristic envelope. On integration 
of Eqs. (11) with respect to z and introduction of symmetry for convenience for 
which u,(x,y,0) = 0, we obtain a two-dimensional "reduced" version of (11), i.e., 


Qyy - S 2 Qxx = S(x,y) 


where 


S(x,y) 


3 

ay 



+ 




- 8 2 



+ 




(13a) 


(13b) 


and the subscript c denotes conditions on the characteristic envelope. Turning 
to the appropriate boundary conditions for (13) , we use (lib) and (12) to obtain 


Qy(X,0) = F(x) 


+ u_ 



-8G'(x) 


(14a) 


where 


S LE 

F(x) = / f. 
0 


xx 


(x,z) dz 


(14b) 


with subscript LE denoting the leading edge value. 

To complete the fonnulation of an initial boundary value problem for Q in 
the x,v plane, we introduce the condition 


Q(-x,|) = o 


(14c) 
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which is appropriate if the parabolic Mach cone intersects the plane of symmetry 
such as for a supersonic leading edge delta wing. Equation (14c) is derived on 
the basis that the cone carries aero lateral load. 

Heuristically, we assert that S=0 in (13b) . It is anticipated that a 
rigorous justification of this assertion can be made using the characteristic 
compatibility relations and the fact that the characteristic surface is an 
equipotential. A similar argument was used to derive the hypersonic integral 
theorems presented in Ref. 5. 

The solution of the initial boundary value problem embodied m (13) and 
(14) is the "simple wave" 

Q = G(x - By) (15) 

where G is obtained from (14) . 


SPECIALIZATION TO DELTA WING 


Referring to figure 14, (15) will be specialized to treat the case m ivhich 
f = 0 in (14b) to illustrate the concepts. Thus, 


t ■ u c — 


on v = 0 


(16) 


and hence 


Q = a(v - x/S) 

where a is to be determined by (16) . 
In particular, 

Q(x,0) = - ax/B . 


60 



CONICAL PROJECTION 
OF MACH WAVE 



Figure 14. Notation for Delta Wing Solution 
by Integral Method 


Hence , 



act 

8 cot A 


where C L is the lift coefficient representing an integral of loads along strips 
such as AA' in figure 13, S^ 2 is the half planform area and A is the sweepback 
angle. Now from figure 14, the Mach wave plane is given by 
z_ = x cot A - y cot A 

and, 

3z c /3y = - cot X 

where X is the ray angle to the point of tangency of the wave to the Mach 
circle. Since u is the two-dimensional value of the "pressure" obtained from 
sweepback theory, evaluation of (16) gives 


jj = sec X cot X = esc X = 2 ctn A 

so that we obtain finally from (17) that 

C L = 4a/ 8 . (18) 

This result checks the appropriate value obtained from the full three-dimensional 
theory. This example illustrates the relative ease that this procedure can be 
applied to obtain an overall force coefficient which avoids entirely function 
theory methods or surface singularity procedures . Extensions to curved leading 
edges and more general geometries appear feasible. 
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SECOND ORDER TRANSVERSE INTEGRAL METHOD 


In this section, we outline the extension of the procedure to treat the 
second order theory problem, which can be written in abbreviated form as 

□u 9 - g(x,y,z) 

u 7 (x,0,z) = h(x,z) 

y 

where g and h are functions of the first order Prandtl Glauert solution. 

Denoting second order quantities corresponding to those previously introduced by 
a subscript 2, and using procedures applied in the previous section, we obtain 


□'Q i = J g dz + S 7 (x,y) 

(19a) 

“ J o 


Q ? (x,0) = - BG ' (x) 
Y 

(19b) 

Q 2 (x,x/B) = 0 

(19c) 


where G 7 is a known function, and 



Introducing the characteristic coordinates 
£ = x - By 
n = x + By , 


Eq. (19a) can be written as 
^2 


3 2 Q- 


3£9n 


= H(5,n) 


( 20 ) 
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and can be solved with simple quadratures using the boundary conditions (19b) 
and (19c). Thus, the complete second order problem for the normal force has 
been reduced to a problem similar to that for the first quantities embodied in 
(13) and (14) , with a slight complication involving inhomogeneities in the strip 
loading equation. In future work, it is intended to solve (19) for a delta 
wing and compare the results for the lift and other forces obtained in this 
fashion with those obtained using the source volume and finite difference 
analyses described previously in Section 4. The obvious advantage of this 
method is that it reduces the three dimensional problem to one in two dimensions 
and has the potential of being generalized to more complicated geometries. On 
the other hand, it has the disadvantage, of yielding-only gross loadings and 
not detailed pressures. 
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6. RESULTS 


6.1 Second Order Theory 

The initial test cases for the source volume formulation of second order potential 
theory were selected to verify the correctness of the theoretical development and 
coding. Available analytic and higher order two dimensional results were used for 
this purpose. Comparisons were subsequently made with measurements for delta wings 
covering the range of leading edge sweep from 50 to 70 degrees, Mach numbers of 2.3 to 
6.0, and a range of angle of attack such that Ma 1. 

Comparisons for a NACA 0012 airfoil at zero angle of attack are presented on 
figures 15 and 16 . The incompressible result of figure 15 established the 
correctness of the second order boundary condition by comparison to the conformal 
mapping solution of Theodorsen. The well known inadequacy of the first order chord 
plane boundary condition transfer in the leading edge region is clearly illustrated. 

The second order result essentially eliminates this deficiency. The compressible 
subcritical comparison of figure 16 was used to establish the correctness of the 
spatial source volume terms by comparison to the higher order Euler equation solution 
of Sells. The results indicate that the source contribution generally improves the 
comparison although the boundary condition correction is the more dominant of the 
two at this Mach number. Higher order contributions are apparently required to 
entirely resolve the leading edge region. 

Numerical second order calculations at M = 4 were performed for a swept taper 
ratio one lifting wing with a 10% circular arc sectinn (figure 17 ) . The results for 
an essentially two dimensional region of the flow were compared with Busemann's or 
Van Dyke's analytic second order expression for yawed wings. 
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The only difference between the analytic and numerical result is in the second 
order net pressure and is due to the chordwise variation of camber induced by the 
source volumes. The camber is measured at the control point, which is at 0.95 chord 
of each vortex panel, while the ACp is assumed to occur at the centroid. For a 
supersonic two-dimensional region the net pressure depends only on the local angle 
of attack. 

The chordwise variation of the three terms which contribute to the second 
order ACp are shown in figure 17c. The magnitude of the three terms is different 
depending upon whether the free stream is placed at angle of attack or the planform 
is place at angle of attack. For either choice in the essentially two dimensional 
region, the exact second order results and the program results reduce to the same 
value. Both show there is no change in lift due to thickness which, is in agreement 
with analytic results for an uncambered two dimensional section which, begins 
and ends with zero thickness. 

Numerical second order calculations for a clipped fifty-five degree leading 
edge sweep delta wing with a 41 circular arc airfoil at M = 2.3 to 4.6 are compared to 
wind tunnel measurements on figures 18 through 20. Predictions of lift surve slope by 
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the first and second order analysis are essentially the same, and both are in 
satisfactory agreement (figure 18) with experiment. The prediction for aerodynamic 
center exhibits a marked improvement for the second order result and is in good 
agreement with test results. A comparison between prediction and measurement as 
a function of angle of attack is presented on figure 19 at a M = 4.6. Figure 20 
presents representative chordwise pressure distribution comparisons . The 
general expression for the upper and lower surface pressures, as predicted by 
second order theory, for a planar wing with zero camber is of the form: 

= b , (~r) + ^ + °<(f) + C j °< 

C r. T Mt) * b t d) - a ,°< " + 

where ap, bp, a.2, b 2 , and Cp depend on planform geometry and thickness distribution. 
The chordwise pressure distributions from the wind tunnel data where curve fit by an 
equation of the form 

C = -f "f ^ ^ ^ + -f 

* o i z. 3 

The data used for the curve fit was limited to ±11° to avoid problems assoicated with 
flow separation. Since there was data available at more than twice as many angles 
of attack as there were unknown coefficients, the curve fit was performed in a 
least square sense. This helped to eliminate the scatter in the wind tunnel data. 

The f Q term represents the Cp due to thickness while the £2 term represents the second 
order contribution to Cp due to lift. The fp term can be used to represent the ACp 
while the fj term was found to be negligible. 
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Calculated second order results for the compression side of a 50 degree swept 
delta wing at a Mach number of 4 and an angle of attack of 5 degrees are compared 
to the higher order method of lines. on figure 21A. Good agreement is indicated and 
a substantial improvement over first order analysis realized. 

Finally, second order results for the compression side of a 70 degree swept 
delta wing at a Mach number of 6 and an angle of attack of 8 degrees are compared to 
measurements on figure 21B. The agreement is satisfactory and again provides a 
substantial improvement over first order predictions. The non-smoothness of the 
second order predictions is attributed to the coarseness of the source volume density 
which for the case under consideration is 10 chordwise by 8 spanwise by 8 vertical. 
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Finite Difference Second Order Theory - Velocity Potential Scalar Approach 

Calculations were performed for flat delta wings at = /2 and sweep angle 
X = arctan (0.7071) with and without fourth-order damping. The results for the 
first order equation are shown in figure 22. Even when points are placed at the 
leading edge, some oscillations are seen inside the Mach cone when no damping 
is used. The addition of damping seems to remove these high frequency oscillations 
as shown in figure 22. Considerable improvements were also noticed when damping 
terms are added for cases with points not on the leading edge of the delta wing. 

Second-order calculations were also performed for flat plate delta wing 
geometry with grid points on the leading ^edge^ Two typical results are shown 
in figures 23 and 24. Considerable improvements over the linear theory 
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Figure 22. Comparison of analytical anil numerical solutions with and without 
fourth-order damping. = / 2 , ?X = arctan 0.7071 




» 3.0 , X = 45° 



(z/x) tan x 

Figure 23. Compression side — delta wing Human boundary 
point operator with damping 



Figure 24. Cocroression side — delta wing Human ooundary 
point operator with damping 
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are obtained when the second-order theory is used. In figure 23, the 
linear and the second-order theory results for H OT = 3, x 3 ^5° are compared 
with Fowell's (Ref. 10) exact solution. At 4° angle of attack, the two- 
dimensional region computed by the second-order theory compares well with 
the exact solution with some deviation inside the Mach cone. In figure 24, 
the exact solution using method of lines (Ref. 11) is shown for comparison 
with the second-order results. For a preliminary design code concept, 
the improvements obtained with the second-order theory over the first-order 
solution is very encouraging at this point. 

One other check was performed to validate the results obtained from the 
computer using second-order calculations. The pressure coefficient for 
a flow turning angle of A8 and Mach number is given by Busemann theory 
(two-dimensional unswept value) 


p - p 

C =* =* ±C. A9 + C_ (A9) 2 ± C, (A9) 3 + ••• 

p K«i 1 3 


(26) 


where A8 is positive when counterclockwise. The upper and lower signs 
refer respectively to left-running and right-running Mach waves, and the 
coefficients are as follows: 


C 


1 


2/VM 2 


1 


C 2 3 [(M 2 - 2) 2 + YM^]/2(M 2 - l) 2 


C 


3 


Mi 


(M 2 - 1) 


7/2 



5 + 7y - 2 y 2 \ 

2(y+1) ) 


3 (Mi - 4/3) 2 
4 (M 2 - 1) //2 ' 

It will be noted that the first term of Eq. (26) represents the linearized 
solution, and the first two terms the second-order solution. 


u. ~ + 28 y 3 + llY 2 - 8 y - 3 

24 (Y+l) 


Equation (26) does not have the sweep effects taken into account. 3y 
replacing M^ by M ot cos x and A8 by arctan (tan A6/cos x) and compensating 
for pressure coefficient normalization with respect to the freestream 
reduced by a cos X factor, we get the sweep effects. For M^ 3 4, 

A8 3 5 a , and x = 50° the first two terms m Eq. (26) (after making the 
sweepback corrections) yield a value of ( c p) secon d order '’-D = 05733831. 

For the same case, the finite differenced second-order calculations give 


(C ) 

P second orde: 


0.0573339. 
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Figure 25 shows a typical wing-body arrangement. Here, a right circular 
cone whose axis is rotated with respect to the wing to avoid a 90° intersection 


at wing-body junction, is placed on top of the delta wing. Figure 26 shows 
another arrangement of wing-body combination with a cosine body with no slope 
discontinuities at the wing-body junction. Calculations performed on these 
wing-body combinations using the velocity potential approach in Cartesian system 
are shown in figures 27 and 28 in terms of constant C^ contour lines. For the 
cone body, the delta wing was at zero degree angle of attack. Thus, in figure 27, 
the Cp in the 2-D region on the wing is zero. For wing-body arrangements once 
the grid points were placed along the leading edge the wing-body junction may 
or may not have points lying on it. In the Cartesian velocity potential 
formulation if the grid points were not 'aligned along the -leading' bdge and' 
wing-body junction the solution tends to be very oscillatory. The spikes that 
are seen in figures 27 and 28 originate at points neighboring the wing-body 
junction and propagate inwards along Mach lines. As expected the cosine body 
produces oscillations of smaller magnitude because of a smooth wing-body junction. 


Finite Difference Second Order Theory — Systems Approach 

Several conical delta wings and wing-body combination cases were run using 
the first and second order system of equations approach. Most of the cases 
used a typical n,£ grid of (11 x 18) and required approximately 400 iterations 
to converge. All the cases were run on the Berkeley CDC 7600 machine.* A 
typical first and second order calculation using the systems approach required 
10 seconds CPU time. 


*This procedure although not strictly speaking necessary was chosen to reduce 
costs during the contractual effort'. 
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Figure 25. ^ Oxc*- « • *”* 
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Figure 


27 Constant Cp Lines for Wing-Rotated Circular Cone Body. 

08 -Y = 65°, a = 0°, Cone Semi-angle = 12.4 , 
tSe Axis of the Cone is^Rotate'd mtfr Respect to Wing 
by 6.2° . 



Figure 28. 


Constant CL Lines on Wing-Cosine Body. - 4, x = 50°, 
a = 5°. F 
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Figure 29 shows the first order solution and -the improvements made using 
the second order solution for a delta wing case of M ot = 4, x = 50° and angle 
of attack = 5° . The system of equation approach yields a much smoother solution 
than the velocity potential formulation, especially for the second order 
calculation. Better flow field resolution in the case of the system of equation 
approach could be due to the conical transformation. 

Figure 30 shows another case of a delta wing for = 6, x = 70° , a = 8° . 

The first and second order solutions are compared with the Gentry/Woodward solution 
presented in Ref. 12, experimental data, the more exact method of lines calculation. 
The second order finite difference solution and the Gentry/Woodward solution 
compare reasonably. It is very clear from this plot that a significant improvement, 
in the results can be achieved using' secorJd order approach. -However some 
disagreement with experiments still remain due to the inadequacy of the linear 
and Gentry theories in predicting the nonlinear shift in conical sonic line 
location. 

Figure 31 shows a typical wing-body calculation for = 6, x = 70°, and 
a = 8°. A conical cosine body having the' following shape (see figure 26) 


= 6x (1 + cos it ■§-) + 


(XX 


0 < 5 < £ b 


" 006 C b < 5 < 5 <ie 

is used, where ^ is the location of body-wing intersection and is the 
leading edge position. The pressure of the body does not change the two- 
dimensional results as long as the body is contained within the Mach cone. The 
boundary condition for the body is also applied at y=0, just like for the wing. 
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Figure 29 . Compression Side — Delta Wing 
4, x - 50°, a - 5° 
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Figure 30. Compression Side — Delta Wing 
M = 6, x - 70° a “ 8° 

OO * /V 



l'lt'.ure ?l. Compression Side - Conical Wing-Boily Geometry 

“u, “ 6, X - 70°, a = 8*< 











Comparing figures 3d and 31 , it can be seen that the influence of the body 
completely changes the second order solution in the three dimensional region. 
Instead of a sonic line in the case of delta wing a compression wave (can be 
a shock also) is formed terminating the two dimensional region. The flow then 
undergoes a gradual compression thereby peaking the C p value. The peak in C p 
occurs around the intersection of the wing-cosine body. On the body the flow 
(cross flow) undergoes a rapid expansion and the C Q value falls down to the 
centerline value. 

Figure 32 shows the results for symmetric subsonic leading edge delta 
wing cases. For the linearized equations there is a square root singularity at 
the subsonic leading edge causing .tcuga-to- infinity . The numerical procedure 

cannot handle this singular behavior properly and hence forms spikes in the 
solution near the subsonic leading edge. Results are shown for = /2, 
semi wedge angle of 3° and for leading edge sweep angles 69° 56' and 51° 21'. 

The first and second order results are also compared with analytical results 
from Ref. 13. The discrepancy between the analytical results and the numerical 
results require further investigation. 
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6.2 Results — Evaluations of Hypersonic Small Disturbance Theory Model 


In this section, we apply the formulation devised in Section 5.1 to treat 
flows over elliptic cones at incidence. This study is intended to illuminate 
accuracy questions in regard to the initial concept and indicate directions 
of refinement prior to application to more general configurations. 

Referring to figure 33, we consider the elliptic cone in the hypersonic 
stream at freestream Mach number and incidence a. In Eq. (7b) in Section 5.1, 
by use of a body axis to wind axis transformation, we find that for the shape 
in figure 33, the body B and shock S have the special representations 


B = r - xf(0) = 0 (la) 

S = r - xg(9) = 0 (lb) 


where r = r/5, y = y/5 = r sin 0, z = 7/5 = r cos 9, r 2 = y 2 + z 2 , and 
9 = tan' 1 y/z. For an elliptical cross section of eccentricity e and major 
axis vertex angle 5, we have, with A s a/5 


f = r/x = 


- A sin 9 + 


V A 2 sin 2 9 + 


s 2 (sin 2 9 


+ e 2 cos 2 


9 ) 


sin 2 9 + e 2 cos 2 e 


(2a) 


e Vsin 2 9 + e 2 cos 2 9 - A sin 9 
sin 2 9 + e 2 cos 2 0 


+ 0(A 2 ) as A -*• 0 


f ' (9) = - (f cos 9) 


f(l - e 2 ) sin 9 + A 
sin 2 9 + £ 2 cos 2 9 + A sin 9 ' 


(2b) 


The pressures on the shock are given by Eq. (3d) where the right hand side is 
calculated using Eqs. (2). Defining this right hand side as a quantity S(9), 
we obtain in the notation of Eqs. (7) and solving (8) 
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f 2 


( 3 ) 



VF 2 + F 2 Vf 2 + £' 2 



For this conical shape, the nonlinear first order partial differential 
equation in two space variables x and 9 becomes the following ordinary differential 
equation for the reduced shock shape g: 

g'(9) = g V(g/T(e)) 2 - 1 . (5) 


Note that Eqs. (3) -(5) apply to an arbitrary cross section conical body, 
a surface integration is used to determine the lift coefficient, then the 
appropriate expression for a conical body .is 

IT 


C T = ff C n- j dS 
L p B 


-I 


c (9) 

TT 

7 


f[f sin 9 - f' cos 9] 

Vf 2 + f ' 2 


d9 


If 


( 6 ) 


In figure 34, results using the present hypersonic small disturbance "wedge" 
model are compared to other theories for the a=0 right circular cone (£=0) case. 

In contrast to the = 00 results previously presented in figure 12, the discrepancies 
increase from 20% at H=0 corresponding to the strong shock = « case to about 
49% at H = 66 = 1 for y=1.4. The discrepancies could be reduced by employing a 
lower value of y associated with a thin shock layer approximation or the previously 
used somewhat inconsistent approximation employed in Ref. 8 of negligible freestream 
pressure compared to static pressure on the body. Further evaluation of this 
case indicates that the present theory makes too gross an approximation of the 
wave angle T by assuming that the magnitude of transverse component of velocity 
|V T | = |q - Ui| is conserved across the shock layer. For a wedge in which the 
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Figure 34. Comparison of Integral Method for 
Normalized Coefficient Pressure 
Over Right Circular Cone at Zero 
Angle of Attack, g = /jz _ 
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flow behind the shock is constant state this assertion is obvious. Accordingly, 
it would be expected that for relatively planar winglike bodies with supersonic 
leading edges this approximation may not produce the serious errors associated 
with the cone. For the latter, however, important flow relief effects occur due 
to the high three dimensionality. Also implicit in the comparison of figure 34 
is the assumption that shock and body pressures are the same. This assertion is 
apparently not as serious as that concerning the wave angle, T. The error in 
the latter is only 20 % but its appearance as a square aggravates the prediction 
sensitivity. 

Prior to discussing directions for refinement, we turn now to the case 
a, e^=0, i.e., an elliptic cone at incidence. With the notation given in 
figure 33, pressures were calculated with our model for an c = 1/2 cone tested 
by Chapkis at M m = S.8 and reported in Ref. 14. Results for the pressures on 
the most lee (9 = m/2) and windward rays (9 = - it/ 2) are shown in figure 35 
for various angles of attack a. For the windward side in which the thin 
constant property shock conditions are most closely met, the agreement is 
reasonable. Degradations occur on the lee side associated with thickening 
of the layer and resulting nonuniformities in the region. In figure 36, the 
experimental meridional variation of pressures is compared for the same cone 
at a - 0 and 2° with the predictions from our model. Serious discrepancies 
occur near the major axis. Newtonian flow results for a * 0 are also shown 
and also indicate serious discrepancies, but near the minor axis. The appropriate 
Newtonian expression for zero incidence is 

2 ■ = c 2 [sin 2 9 + e 2 cos 2 8]/{[sin 2 9 + c 2 cos 2 8] 

(1 - e 2 ) sin 2 cos 2 8} . (7) 
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Figure 36. Comparison of Approximate Shock Cp with Experimental 
Body C - Elliptic ('one 



For the prediction of by a surface pressure integration, one might 
expect a beneficial error reduction due to the tilting of the pressure vector in 
the region where it contributes the least to the lift and the error is the 
greatest. Unfortunately, for an ellipse, this tilting is also localized to the 
major axis and no substantial accuracy improvement occurs, as shorn m the 
meridional loading curve of figure 37. Using Eq. (6), we obtain a of 0.05056 
in contrast to a value of 0.0278 using Chapkis' measured pressures. 

As indicated previously, a crucial feature of the model is its capability 
to accurately predict the shock wave shape. In the foregoing phase involving 
the determination of C^, the application of (5b) in Section 5.1 to the 
determination of was considered on the grounds that compensating errors 
could be involved in the evaluation of the crv product and the estimation of LA. 
The determination of AA involves the solution of (5) for g. Cauchy-Euler 
marching was attempted, but it became evident that the singular behavior near 
0 = - j will have to be better understood before this scheme is successful. 

More accurate prediction of the shock shape for this and generalized 

geometries will require suitable refinements in the model in which presumably, 

more localized integral forms of the equations than (5) in Section 5.1 are 

used. The existing model utilizes onlv the surface boundary conditions but 

does not incorporate the differential equations, particular lv that of continuity 

This aspect is particularly important in distinguishing the highly three- 

dimensional conical flow considered previously within the framework m the 

quasi- two- dimensional assumption of wJ = V T . In future studies, finite 

1 ! S 1 B 

volume procedures involving these conservation laws will be more fully developed. 
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7. CONCLUSIONS 


Based on the theoretical development and comparison with higher order 
results/experimental measurements, the following findings and conclusions are 
made. 

1. Second order potential theory provides a systematic means of 
extending linear theory to values of the similarity parameter 
M6 of order one. 

2. Improved prediction of supersonic/hypersonic aerodynamic characteristic 
and surface pressures for simple three dimensional shapes has been 
demonstrated for numerical second order analysis. 

3. Additional effort is required to extend the non-linear potential analysis 
to more general wing-body arrangements and refine the treatment of 
subsonic edges. 

4. Exploratory effort utilising hypersonic small disturbance theory 
integral techniques indicates that the approach requires additional 
development prior to obtaining a rapid gross aerodynamic prediction 
method. 

5. The systems approach has a strong potential for any future code 
development to handle general three dimensional configurations in much 
the same manner as existing Euler codes but with sufficient computer time 
and cost savings to permit use as a preliminary design tool. 
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